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ABSTRACT 

The discovery of extra-solar planetary systems with multiple planets in highly eccen- 
tric orbits (~ 0.1 — 0.6), in contrast with our own Solar System, makes classical secular 
perturbation analysis very limited. In this paper, we use a semi-numerical approach 
to study the secular behavior of a system composed of a central star and two massive 
planets in co-planar orbits. We show that the secular dynamics of this system can be 
described using only two parameters, the ratios of the semimajor axes and the planetary 
masses. The main dynamical features of the system are presented in geometrical pic- 
tures that allows us to investigate a large domain of the phase space of this three-body 
problem without time-expensive numerical integrations of the equations of motion, and 
without any restriction on the magnitude of the planetary eccentricities. The topology 
of the phase space is also investigated in detail by means of spectral map techniques, 
which allow us to detect the separatrix of a nonlinear secular apsidal resonance. Finally, 
the qualitative study is supplemented by direct numerical integrations. Three different 
regimes of secular motion with respect to the secular angle Azu are possible: they are 
circulation, oscillation (around and 180°), and high eccentricity libration in a nonlin- 
ear secular resonance. The first two regimes are a continuous extension of the classical 
linear secular perturbation theory; the last is a new feature, hitherto unknown, in the 
secular dynamics of the three-body problem. 

We apply the analysis to the case of the two outer planets in the v Andromedae 
system, and obtain its periodic and ordinary orbits, the general structure of its secular 
phase space, and the boundaries of its secular stability; we find that this system is 
secularly stable over a large domain of eccentricities. Applying this analysis to a wide 
range of planetary mass and semimajor axis ratios (centered about the v Andromedae 
parameters), we find that apsidal oscillation dominates the secular phase space of the 
three-body coplanar system, and that the nonlinear secular resonance is also a common 
feature. 
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1. Introduction 

At the present time, at least ten sun-like stars in our galactic neighborhood are known to 
harbor planetary systems with two or more planets^. Two (perhaps three) of these include two 
planets locked or close to orbital mean motion resonances, and in these cases the orbits of the 
planets are significantly perturbed on orbital timcscalcs by the mutual gravitational interactions 
of the planets. In all other non-resonant systems, the long term orbital variations are qualitatively 
and quantitatively described primarily by the secular effects of the gravitational interactions of the 
planets. The secular interactions can produce interesting orbital dynamical effects, such as align- 
ment (or anti-alignment) of the lines of apsides and large amplitude variation of the eccentricities. 
The planetary system of v Andromedae has received particular attention in this regard, owing to 
the very small difference in the arguments of periastron and the surprisingly large eccentricities 
of its two outer planets (Chiang et al. 2001, Malhotra 2002, Chiang and Murray 2002). An un- 
derstanding of the secular dynamics may potentially help identify the origin of the large orbital 
eccentricities of some extra-solar planets (Malhotra 2002). 

For the planets of our solar system (excluding Pluto), there exists a well-developed theoretical 
framework, based upon the classical Laplace-Lagrange perturbation theory (Laplace 1799; see, e.g., 
Murray and Dermott 1999) for the analysis of the secular orbital perturbations of the planets. 
However, in contrast with our solar system, the known exo-planetary systems have generally larger 
planet masses, smaller semimajor axes, and larger orbital eccentricities; these factors, particularly 
the larger eccentricities, limit the application of classical linear secular perturbation theory for these 
systems. 

In the present paper, we analyze the secular dynamics of the three body problem consisting of 
a central star and two coplanar massive planets with large orbital eccentricities. Our approach is, in 
essence, the nonlinear generalization -to arbitrary eccentricities- of the Laplace-Lagrange secular 
perturbation theory. We use a semi-numerical approach, which employs a numerical averaging of 
the short period gravitational interaction of the planets, to determine the interaction Hamiltonian 
describing the secular perturbations in the system. We then use this secular Hamiltonian and a 
Hamilton-Jacobi approach to discern the geometry of the phase space and the main dynamical 
features of the system. We show that the secular phase space structure of the three body (two 
planet) system is determined by only two parameters, the ratio of the planet masses and the ratio 
of their orbital semimajor axes. We calculate the periodic orbits, their stability characteristics, and 
the domain of oscillation about the periodic orbits. We present our analysis over a wide range of 
planetary mass and semimajor axis ratios, including those pertinent to the v Andromedae planetary 
system. 

In the low-to-moderate eccentricity regime, there exist two periodic orbits which can be iden- 
tified with the aligned and anti-aligned modes of linear secular theory (the alignment refers to the 



^http://www. obspm.fr/planets 
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longitudes of periapse of the two planets). Small amplitude oscillations about these periodic orbits 

are sometimes referred to as 'secular resonance' or 'apsidal resonance' in the literature. However, 
in the present work, we reserve the "secular resonance" terminology for a different type of motion, 
for the following reason. In the regime of moderate-to-high eccentricities, we find that the aligned 
mode bifurcates into a pair of stable and unstable periodic orbits, and a new, third, regime of 
motion becomes possible. To the best of our knowledge, this is a previously unknown feature of 
the secular dynamics of two-planet systems. It consists of a zero-frequency separatrix bounding a 
nonlinear resonance zone in the secular phase space at large values of the eccentricities. [An anal- 
ogous feature in the secular dynamics of massless particles in the solar system has been described 
previously in Malhotra (1998).] We verify and supplement our analysis of the numerically averaged 
system with direct numerical integrations of the unaveraged system. 

The paper is organized as follows. In section 2, we describe our model, the numerical averaging 
approach, and the geometrical method for locating the periodic orbits of the system using the global 
constants of motion. In section 3, we apply this approach to obtain a picture of the secular phase 
space neighborhood of the outer two planets of v Andromedae. In section 4, we explore the secular 
phase space structure over a wider range of the parameter space of planetary masses and semimajor 
axis ratios. Finally, we provide a summary and discussion of our results in section 5. 



2. Secular dynamics in the three-body problem 
2.1. The model 

Consider two planets of masses mi and 1712 orbiting a central star of mass M. Hereafter, the 
index i = 1, 2 stands for the inner and outer planet, respectively. In the heliocentric reference 
frame, the set of canonical variables introduced by Poincare (1897) consists of the planets' position 
vectors fi relative to the star and their conjugate momenta pi = rrii where pi are the position 
vectors relative to the center of gravity of the three-body system (Yuasa and Hori 1979, Laskar 
and Robutel 1995). In these coordinates, the Hamiltonian H of the three-body problem may be 
written in the form 

Pi Mi^N G mi 1712 , (Pi-P2) 

direct part indirect part 

Keplerian part 

where G is the gravitational constant, Hi = G{M + mj), m[ = rrii M/{M + rrii) and A = |ri — r2|- 
The first term produces Keplerian motions of the planets; the second and third terms produce 
direct and indirect perturbations among the planets, respectively. 

Associated with the Keplerian part of the Hamiltonian, a set of mass-weighted Poincare elliptic 
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variables is introduced as 

Ai=mean longitude, Lj = m[.yj]liai 



—■ci7j=longitude of perihelion, Li — Gi= — -y/l — e|), 

— rij=longitude of node, Gi — Hi=Li\ 1 — e? (1 — cos ii 



where Oj, Cj and are the heliocentric canonical semi-major axes, the eccentricities and the inclina- 
tions of the planets, respectively. For the v Andromedae system, dynamical stability considerations 
suggest mutual inclinations less than 20° (Stepinski et al. 2000, Chiang et al 2001). Hence, it is 
reasonable, in first approximation, to set the planetary inclinations to zero. The Hamiltonian (1) 
can then be written, in terms of canonical elliptic variables, as 

^ = - E ^yf - X ^(^- - ^^)' (2) 

1=1 ' 

where the sum describes Keplerian motions of the planets and R is the disturbing function. 

To study the secular behavior of the system, we perform a canonical transformation to the set 
of the new variables: 

Al Li 

(3) 

Acij = CC72 — vji K\ = Li — G\ 

-W2 K2 = {Li - Gi) + (L2 - G2). 

Then we apply the following averaging procedure to the Hamiltonian system of Eq. (2). The 
averaging is done with respect to the mean longitudes of the planets A, and removes from the 

problem all the short periodic oscillations. This averaging is done by a numerical process through 
the double integration of the function (2). Because the Keplerian part does not depend explicitly 
on the angles or on Ki, it contributes only an inessential constant to the secular Hamiltonian, and 
can be neglected; thus, the secular Hamiltonian is defined by 

^sec =-7:^2 r r X ^d^^d^-- 

In practice the numerical integration need be done over only the direct part of the disturbing 
function, because the indirect part, Eq. (1), does not contain secular terms (Murray &; Dermott 
1999). The value of the Hamiltonian and its derivatives are computed using Eq. (4) for any 
given point of the phase space. We emphasize that, in the present work, we do not perform the 
developments of H in the power series of the planetary eccentricities nor in Fourier series of the 
angular variables, such as in many previous analyses (e.g., Brouwer and Clemence 1961). This 
means that the secular motion of the planets is described very precisely with our method, without 
any restriction about the magnitude of their eccentricities. The conditions for the applicability of 
the analysis are that the system is sufHciently far from a strong mean motion resonance and that 
the orbits are assumed to be coplanar. 
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After the elimination of the short periodic terms, the averaged Hamiltonian ifgec does not 
depend on Ai and A2; consequently, Li and L2 (thus ai and 02) are constant in time and serve 
simply as parameters in i^scc- In addition, due to the D'Alembert's rule, the tu-dependcncc in the 
averaged disturbing function is only through Aw. So, —W2 is a cyclic coordinate in ifseo hence 
its conjugate momentum K2 is also a constant of motion. Thus, the averaged planar three-body 
problem is separable and integrable. Hamilton's equations for this problem are as follows: 

and 

K2=0; ^2=-fc. (6) 

The solution for K2 is trivial, K2 = -fC| = const. This means that the averaged planar three-body 
problem has been reduced to a one-degree-of-freedom dynamical system (5), in which K2 (and Li 
and L2) are constant parameters. One important conclusion that follows is that the variations of 
Aw and Ki are described by a single frequency. ^ The solution for Ki and Aw can be obtained 
by integrating simultaneously the equations (5), for given fixed K2 (and Li and L2). Next, the 
equation (6) of motion for W2 can be integrated through a simple quadrature, using the fixed 
value of K2 and the obtained functions Ki{t) and Aw{t). Note that this reveals the second proper 
frequency in the planar secular system. Thus, the general solution of the full two-degrees-of- freedom 
dynamical system described by i^scc in Eq. 4, is quasi-periodic, with two proper frequencies. We 
will see that particular solutions, when the two planets precess in lockstep with fixed eccentricities 
and aligned (or anti-aligned) apsidal lines, are periodic solutions of the full problem, but are fixed 
points of the reduced one-degree-of-freedom system, Eq. (5). 

After obtaining the solution in terms of the canonical variables, we can use the expressions in 
Eq. (3) to do the inverse transformation to the orbital elements to describe the secular time vari- 
ations in the individual eccentricities and apsidal longitudes of the two planets. To represent the 

results of our study in familiar orbital elements, rather than the canonical set (Ki, Aw; K2, —W2), 
we will use the particular set of variables, ciexpiAw and e2expiAti7. This set of variables is 
preferable over the classical variables (ei expiwi and e2expiro2) as the qualitative aspects (align- 
ment of lack thereof of the apsides) are more readily obtained with these variables. To obtain the 
time variations of these, we need to solve only the first set of Hamilton's equations, Eq. 5. The 
possible solutions are either fixed points or periodic solutions, characterized by a single frequency. 
However, as noted above, solving the equation of motion for W2 (Eq. 6) by quadrature and the 
transformation to the orbital elements, allows us to return to the full two-degrees of freedom system 
and the classical variables eiexpiroi and e2expitU2. This full system has two proper frequencies. 
Keeping this in mind, we refer hereafter to the fixed points of Eq. 5 as periodic solutions (of the 
full problem), and to the general solutions of Eq. 5 as quasi-periodic (or ordinary) solutions. 



^This result can also be derived from the classical Laplace-Lagrange secular theory for two planets. 



(5) 
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The classical Laplace-Lagrange planetary theory (see Murray and Dermott 1999) describes 
well the behavior of the secular system, for small values of the planet eccentricities. The approach 
used in this work is an extension of the linear approximation to the domain of the moderate-to- 
high eccentricities. Hence it is expected that the main features of secular motion provided by the 
classical theory will be observed. Over the whole paper, we will compare the results obtained to 
those derived from the Lagrange-Laplace approach, in order to determine the domains of validity 
of the linear approximation. 

The main features of secular motion of planetary systems are defined by such global quan- 
tities as total energy (E), total angular momentum (TAM), and Angular Momentum Deficit 
(AML»)(Poincare 1897, Laskar 2000): 

TAM = Li + L2-k2 
AMD = K2. 

Due to the fact that, in the secular problem, semimajor axes are constants of motion, AMD and 
TAM are equivalent. AMD has a minimum value (zero) for circular orbits and increases with 
increasing eccentricities. The behavior of TAM is reverse of that of AMD. 

In order to apply the developed model to the planets C and D of the v Andromcdac system, we 
adopt the following parameters: the masses M = l.SMgun, rrii = 1.83Afjup and 1112 = 3.97Mjup, for 
the central star, inner planet C and outer planet D, respectively; the semi-major axes ai = 0.83 AU 
and a2 = 2.5 AU. These assume that both planetary orbits are observed edge-on, so sini = 1 for 
each planet, where i is the orbital inclination to the plane of the sky. Although the observational 
data do not constrain the orbital inclinations, dynamical arguments suggest that the uncertainty 
in the individual planetary masses is less than a factor of ~ 2 (Stepinski et al. 2000, Chiang et 
al. 2001). As we will see in the next section, of all these parameters, only two are important 
for the phase space structure: the semimajor axes ratio, ai/a2, and the mass ratio, mi/m2- The 
description of the secular behavior of this system is completed by the initial values of ei, 62 and Am. 
We adopt the following updated values (quoted in Chiang & Murray 2002): ei = 0.252, 62 = 0.308 
at Aw = 0; formal uncertainties in the eccentricities are on the order of 20%. However, we must 
be mindful that the observational uncertainties in orbital parameters of the known multiple-planet 
extra-solar systems are not well determined (Ford 2003). 



2.2. Energy and AMD Levels on the Representative plane 

In this section we present a qualitative geometrical analysis of the Hamiltonian system, Eq. (4). 
We begin this analysis by plotting the energy and AMD level curves in the space of initial conditions. 
The space of initial conditions of the Hamiltonian given by Eq. (4) is three-dimensional, but the 
problem can be reduced to the systematic study of a representative 2-D plane instead of the 
3-D space as follows. In the integrable secular problem, the angular variable Aw can either 
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circulate or oscillate (about or 180°). In both cases, it goes through either or 180° for all initial 

conditions. Hence, without loss of generality, the angular variable Atu can initially be fixed at one 
of these values. The space of initial conditions can then be presented in the (ei,e2)-plane of initial 
eccentricities, where the initial value of the angle Aro is fixed at either zero or 180°. 

Figure 1 shows the levels curves of the energy and AMD (hereafter denoted as i^2^1evcls) on 
the (61,62) representative plane. These level curves were constructed using the parameters of the 
V Andromedae system of planets C and D (assuming edge-on coplanar orbits). The energy levels 
are shown by solid lines, whereas the if2-levels are shown by dashed lines. The domain above the 
thick line at high eccentricities is the region where the planetary orbits cross each other; since there 
are no mechanisms which can protect the planets from close approaches between them (in contrast 
with the case of mean motion resonances), this domain was excluded from our studies. Finally, the 
location of the v Andromedae system (planets C and D) is shown by a full circle, with uncertainties 
of 20% for the planetary eccentricities. 

The structure of the energy levels is independent of the planetary masses; it depends upon 
only the ratio of the semi-major axes, ai/a2- Indeed, as can be seen in Eq. (4), the secular energy is 
directly proportional to the product of the planetary masses and the gravitational constant. Since 
the term 1/A in the direct part of the disturbing function R is independent of planetary masses 
and the indirect part does not contain secular terms, the secular energy can be normalized by the 
factor ^il^L™2._ Then the normalized secular energy, and thereby the structure of the energy levels, 
does not depend on the individual values of the semimajor axes, but depends only on the ratio, 
ai/a2- 

The geometry of the i^2^1cvcls, however, depends on both the semimajor axes and the masses. 
Prom the expression for K2 given by Eq. (3), we can see that a natural normalization factor for 
this parameter is the constant L2. This normalization reduces the semimajor axes dependence to 
only on the ratio, 01/02. The mass dependence of K2 goes as ~ ^(S^)^/^- Therefore, for 
rrii <^ M, in the first order mass approximation, the mass dependence of K2 can be reduced to 
only on the mass ratio of the planets, mi/ni2- For example, in the v Andromedae system, the 
planetary masses are of order 10~^ (in units of the stellar mass), and the first-order approximation 
works sufficiently well in this case. On the other hand, in the case of a binary stellar system with 
one planet, for example, where the mass of the companion is comparable to the mass of the central 
star, this approximation is not valid. For the calculations presented in this paper, we do not make 
this approximation, but rather adopt the values of M, mi, m2, oi, 02 of the v Andromedae system 
stated in the last paragraph of section 2.1. 

2.3. Geometrical pictures of the secular behavior 

In this section, we proceed with a qualitative study of the developed model and give a simple 
geometrical interpretation of the main features of secular motion. It should be emphasized that all 
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Fig. 1. — Energy and K2 - levels of the secular Hamiltonian given by Eq. (4) on the (61,62) 
representative plane of initial conditions. The energy decreases whereas K2 increases with increasing 
|6icosAt:i7| and 62. The signs + or -, preceding the variable eiCOsAtu, indicate that the initial 
values of Azu are zero or 180°, respectively. In constructing the level curves, we adopted the 
following values for the parameters: planet mass ratio 1111/1712 = 0.483 and semimajor axes ratio 
ai/a2 = 0.332 (corresponding to those of the v Andromedae planets C and D). The empty region 
above the thick line is the domain of crossing planetary orbits. The initial conditions of the v 
Andromedae system in this plane is indicated by a full circle symbol; nominal 20% uncertainties 
are indicated for the eccentricities by error-bars. 
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information available from this study is obtained without integration of the equations of motion, 
Eq. (5). The study is based on the fact that the secular energy Hsec and AMD (i.e. K2) are both 
conserved along one secular path. This is reflected on the representative plane of initial conditions 
in the following way: Independent of whether the motion is circulation or oscillation of Aw, an 
orbit generally passes twice through Aw = or 180°. (The exceptional case of singular orbits and 
periodic orbits will be described below.) This means that an individual orbit is represented by two 
points on the (ei,e2) plane: one point being the chosen initial condition, and the other point being 
its counterpart. Due to the energy and AMD conservation, the counterpart of one initial condition 
can be found easily as an intersection of the corresponding energy and iC2-levels. 

Figure 2 shows two representative (ei,e2)-planes of initial conditions which are identical to 
the plane shown in Fig. 1, except that only one energy level is plotted on each of them by thick 
curves. In the top panel, the energy level corresponds to the normalized secular energy of the 
"real" v Andromedae system {Hsec = —0.41391), while, in the bottom panel, it corresponds to a 
lower value Hs^c = —0.4295. We will first concentrate our attention on the top panel. The location 
of the V Andromedae system on the (ei,e2)-plane is shown by a full circle (20% uncertainties are 
indicated by error-bars). The ii'2-lcvel corresponding to the system is shown by a dashed curve. 
This i^2^evel crosses the energy level at two points: one point is the initial condition marked by 
the full circle symbol, and other point is a counterpart of the first point. The two points belong to 
the same solution of the Hamiltonian system described by Eq. (4), and represent a quasi-periodic 
variation in the eccentricities and in Aw. From the condition K2 = const, we can see that the 
secular oscillations of the eccentricities always occur with opposite phases: when the eccentricity 
of the outer planet 62 is minimal, the eccentricity of the inner planet ei is maximal. Conversely, 
when the eccentricity 62 is at a maximum, ei is at its minimum. 

The qualitative behavior of the angular variable Aw also can be obtained from the geometrical 
analysis of the representative (ei,e2)-plane in Figure 2 top. Indeed, if both points of an orbital path 
are located at the positive right-hand side of the (ei,e2)-plane, the angle Aw oscillates around zero. 
Conversely, if both initial conditions are located at the negative left-hand side of the plane. Aw 
oscillates around 180°. Finally, when two initial conditions are at the opposite half-planes, the 
corresponding orbits of both planet are circulating with respect of the angle Aw. The eccentricity 
of the more massive outer planet 62 is minimal (consequently, ei is maximal), when Aw = 0, and 
is maximal (ei is minimal), when Aw = 180°. In this way, from the geometrical picture presented 
in Figure 2 top, we can calculate the ranges of the eccentricity variation for both planets of the 
system under study. 

The analysis of the geometry of the energy and ii'2-levels on the (ei,e2)-plane described above 
is applied for any generic initial condition. However, for a given energy level, there are four 
characteristic K2-levels which define the phase space structure and the qualitative character of 
the secular solutions. These are shown in Figure 2 top by continuous curves. The topmost K2- 
level is tangent to the energy level at the point marked as P"*". The lowermost i^2^1cvel has a 
minimum possible value and is tangent to the energy level at the point marked as P~ . These two 
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Fig. 2. — Top: Energy level with Hsec = —0.41391 (thick curve) and four characteristic i^2^cvcls 
(thin curves) . The location of the current v Andromedae system is indicated by a full circle symbol 
with 20% uncertainties error-bars. The points P"*" and P~ show the location of oscillation centers, 
around zero and 180°, respectively. The initial conditions that correspond to oscillation around 

zero are situated on the energy segment between the points A and B. Those with oscillation around 
180° between the points C and . Bottom: Energy level for Hgcc = —0-4295, and the same four 
characteristic K2 levels as in the top panel, plus an additional K2-level which shows the location 
of a new libration center E and its two counterparts, F and G. Note the different scales of the axes 
in the two panels. 
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K2 levels represent extreme values of K2 for the given energy level. The orbital solutions associated 

with these two extrema correspond to periodic solutions of the secular Hamiltonian system. These 
periodic orbits are characterized by constant values of the eccentricity and of Atu. 

The other two characteristic i^2^1evels contain the initial conditions ei = and 63 = 0. One 
of these passes through the points A and B in Fig. 2 top and other one passes through the points 
C, and D~ . Note that, at 62 = the angular variable Atz7 is undefined. The points and 
D~ both represent the initial condition 62 = on the (61,62) representative plane. 

All initial conditions along a given energy level on the representative plane can be classified 
according to the motion of the corresponding angle Aro. The periodic solution at point is 
characterized by zero amplitude oscillation. The orbits which start on the segment AB of the 
energy level are characterized by oscillation of Aro around zero. The center of oscillation is the 
point P"*". The range of the eccentricity variation is given by the projections of the initial conditions 
on the corresponding axes. The periodic solution with zero amplitude oscillation around 180° is 
at P~ . The orbits which start on the interval between CD^ are oscillating around 180°. Once 
again, the planetary eccentricities vary in the range given by the projections of initial conditions 
on the corresponding axes. Finally, initial conditions in the intervals CA and BD^ correspond to 
circulating orbits with respect to Aro. 

Hereafter we designate the solutions near P+ as Mode I and those near P~ as Mode II. 
Near Mode I, the planetary apsidal lines are aligned, the secular angle Acc7 oscillates about 
and the planet eccentricities undergo small oscillations about the values of the Mode I periodic 
solution. Near Mode II, the apsidal lines are anti-aligned, Atu oscillates about 180° and the planet 
eccentricities oscillate about the values of the Mode II periodic solution. Thus, the two modes 
describe two opposite stable ways of the planetary system to be aligned. Between the domains of 
Mode I and Mode II, there is a region of the phase space in which the motion of the angle Aro is a 
circulation. The secular oscillations of the planet eccentricities always occur with opposite phases. 
The eccentricities undergo secular variations in such way that, when Atu = 0, the eccentricity of 
the outer planet is minimal and one of the inner planet is maximal. Conversely, when Atu = 180°, 
the eccentricity of the outer planet D is maximal and of the inner planet is minimal. 

The bottom panel in Fig. 2 presents by the bold line the level curve of a lower secular energy. 
In this case, the domain covers higher eccentricities, and even eccentricities close to 1. All points 
P"*", P^, A^ B, C, and have the same interpretation as discussed above for the top panel. 
In addition, this figure exhibits a new feature which is the existence of a fifth characteristic K2- 
level shown by points E, F, and G. The new ii'2-level is tangent to the energy level at point E, 
and crosses it at points F and G. This situation means that there is one orbit which crosses the 
61 cos Azu axis at three different points. To understand the behavior in this complicated case, we 
will need to integrate the equations of motion, Eq. (5); the results of the integrations are presented 
in the next section. 



- 12 - 



3. Secular Dynamics of the v Andromedae system: Periodic and Ordinary 

Solutions; Nonlinear Secular Resonance 

In this section, the insight gained by examining the geometrical structure of the levels on the 
representative plane is presented in quantitative form for the system of v Andromedae planets C 
and D. Figure 3 shows the locations of periodic orbits and the various domains of oscillation or 
circulation of Azu in the (ei,e2)-plane of initial conditions, constructed for the parameters of this 
system. The detailed description and interpretation of this figure is as follows. 

The upper boundary of the domain where the planetary orbits do not cross each other is shown 
by the thick line on the (ei,e2)-plane in Fig. 3. This curve is obtained from the simple geometrical 
condition of two crossing orbits, such as: 62 = 1 — ai(l — ei)/a2 and 62 = ai(l + ei)/a2 — 1, for 
Aw = 0, and 62 = 1 — ai(l + ei)/a2, for Aw = 180°. Two families of stable periodic orbits are 
shown by the the thin continuous curves: one is located in the positive half-plane {Aw = 0), and 
other in the negative half-plane {Aw = 180°). The two half-planes are separated by the dashed 
vertical line representing the initial condition ei = 0. 

The periodic solutions of the Hamiltonian Hgec given by Eq. (4) are defined by conditions: 

The Hamiltonian iJsec is an even function of Aw. Hence, the first of these equations provides 
immediately the trivial solution Aw = (mod tt), for non- vanishing Ki and Ki 7^ K2. The second 
equation from Eq. (7) is solved numerically for the two cases. Aw = and Aw = tt, and the 
condition K2 = const. The stability of the obtained periodic solutions is defined by the behavior 
of the Hessian matrix of -ffsec evaluated at the periodic orbit. 

Recall that in linear secular theory, the general solutions are a linear superposition of two linear 
eigenmodes (e.g., Murray Sz Dermott, 1999). Our semi-numerical analysis yields the nonlinear 

generalization of the linear secular theory. In the positive half-plane of initial conditions in Fig. 3, 
the domain of Mode I, i.e. oscillation around zero, is defined by all initial conditions between the 
vertical line ei = and the dashed curve. The domain below the dashed curve is of circulating 
orbits. The domain of Mode II is in the negative half-plane: the initial conditions for orbits 
oscillating around 180° lie between the horizontal axis 62 = and the dashed curve. The domain to 
the right of the dashed curve and up to ei = is of circulating orbits. In the limit of linear secular 
theory, the two modes would be represented in Fig. 3 by straight lines tangent at the origin to the 
continuous curves that are the nonlinear Modes I and II. It is clear that the departures from the 
linear secular theory become significant for ei > 0.3. It is interesting to note that the nonlinearity 
is not symmetric for the two modes; it is especially severe for Mode II even at relatively small 
values of the eccentricity, 62 ^ 0.1. 

Finally, the nonlinear analysis of the secular Hamiltonian yields a new feature not previously 
known, namely unstable periodic orbits at moderate-to-high eccentricities. In Fig. 3, these cor- 
respond to the boundaries of the shaded region near the upper right corner. These orbits were 
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Fig. 3. — Periodic orbits and domains of oscillation and circulation in the (ei, 62)— plane of initial 
conditions, obtained for the parameters of the v Andromedae system. The continuous thin curves 
plot the location of the stable periodic orbits: Mode I in the positive half-plane, Mode II in 
the negative half-plane. The domains of oscillation of Aw about (Mode I) and about 180° 
(Mode II) are bounded by the dashed lines; the regions of circulation of Am are also labeled. The 
shaded region at upper right is a family of unstable periodic orbits. The location of the current v 
Andromedae system (planets C and D) is shown by the full circle symbol. 
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previously identified in Fig. 2 by the points E, F, and G. Their interpretation requires additional 
information, which we obtain by integration of the equations of motion Eq. (5). We describe this in 
more detail below, where we show that the shaded region in the (ei,e2)-plane represents a zone of 
nonlinear secular resonance, consisting of librating orbits bounded by an infinite-period separatrix. 

The location of the current v Andromedae system (planets C and D) in the ei,e2 plane of 
initial conditions is shown by the full circle symbol in Fig. 3. For the presently known best-fit 
planetary parameters, this system is within the Mode I domain with Aw oscillating about zero, as 
has been noted in previous studies (Chiang et al. 2001, Malhotra 2002, Chiang & Murray 2002). 

Figure 4 displays the solutions of the equations of secular motion, Eq. (5), in the two eccentricity 
planes, (cj cos Att7, sin Aro). The top panel shows the secular variations of the inner planet and 
the bottom panel shows the secular variations of the outer planet. These solutions were obtained 
by numerical integration of a large set of initial conditions along the energy level corresponding 
to the V Andromedae system (sec top panel in Fig. 2). In each plane we see two fixed points, 
corresponding to the Mode I and Mode II periodic solutions. The ei plane is dominated by Mode I, 
whose fixed point lies in the central region; the Mode II fixed point lies near the left-hand boundary 
of the energy manifold. Conversely, in the 62 plane, the Mode II periodic solution lies in the central 
region, while the Mode I periodic solution lies near the right-hand boundary of the energy manifold. 

The smooth curves surrounding each of the two fixed points are periodic solutions of the re- 
duced one-degree-of-freedom system, but quasi-periodic solutions of the full two-degrees-of-freedom 
secular system; as noted previously, we refer to these as simply quasi-periodic solutions. Note that 
even though the motion of the angle Aru may be either oscillation (about or 180°) or circulation, 
there is no separatrix associated with an unstable infinite-period solution. To better understand 
this feature, we plot by red and blue curves two particular solutions in Fig. 4. These solutions 
are associated with the singularities in Eqs. (5), which occur at Ki = and Ki = K2 or, in the 
eccentricity notation, at ei = and 62 = 0. (As is well known, these are not physical singularities, 
but merely due to the choice of variables.) The solution shown by the red curve was obtained with 
initial condition ei ~ and is seen as a smooth curve passing through the origin in the ei plane 
(upper panel). In comparison, the solution obtained with initial condition 62 ~ and shown by 
the blue curve, appears as the 'false' separatrix in the ei plane, separating the domains of different 
modes of motion. An analogous situation is seen in the 62 plane (lower panel), where the blue 
curve is a smooth curve passing through the origin and the red curve separates the domains of 
two modes of motion. We emphasize that the boundary separating the oscillation of Avj about 
either or 180° is not a zero-frequency separatrix, but simply the quasiperiodic solutions that pass 
through the singular initial values, ei = or 62 = 0. The separatrix-like curves in Fig. 4 are owed 
to the effects of projecting the canonical solutions on to a plane of non-canonical variables; these 
are better visualized over a sphere (for more details, see Pauwels 1983). 

Finally, the solution corresponding to the current best-fit initial conditions of the v An- 
dromedae planets C and D is shown by the green trajectories in Fig. 4. This solution oscillates 
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Fig. 4. — Secular phase space of the v Andromedae planets C and D, at -ffgec = —0.41391. Locations 
of the Mode I and Mode II are given by the fixed points on the top and bottom panels. The red 
and blue curves represent orbits passing through the ei = and 62 = singularities of Eq. (5). The 
secular variations of the v Andromedae planets C and D for the present best-fit initial conditions 
are shown by green curves (see text). 
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Fig. 5. — Spectral map corresponding to the solutions shown on the planet C panel in Figure 4. 
The abscissas (ei cos Aw) in the spectral map are the initial conditions at Aw = and Aw = 180°. 
The dashed line indicates the initial value of ei cos Aw of the v Andromedae system. 



about the Mode I periodic solution; the oscillation amplitude of Aw is 25° about 0. Note, however, 
that the oscillation amplitude is quite sensitive to the initial conditions. This is most readily seen 
in the structure near the green curves in Fig. 4: small changes in initial conditions lead to large 
changes in the oscillation amplitude. A complete analysis of the observational uncertainties in 
the parameters of this system is beyond the scope of the present work. Suffice to note that these 
uncertainties (see Fig. 3) do not exclude significantly larger or significantly smaller amplitudes of 
oscillation of Aw. For example, assuming all other parameters are fixed, the 20% nominal uncer- 
tainty in the eccentricities implies that the Aw oscillation amplitude may be as large as 70° or as 
small as 9°. As we discuss in section 5, this uncertainty is problematic for theories of the origin of 
the large eccentricities. 

Figure 5 shows the spectral map obtained from the Fourier analysis of solutions calculated 
with initial conditions on the horizontal axis on the planet C panel in Fig. 4 top. The map shows 
the frequency of all peaks present in the Fourier spectra of the eccentricities ei and 62 that are 
clearly above noise level (arbitrarily defined as 20 per cent of the largest spectral peak). Besides the 
discussion above about the nature of the solutions separating the oscillation zones, one may note 
the continuous evolution of the frequencies on the corresponding spectral map. The map shows only 
one frequency and its higher harmonics that is characteristic of one-dcgrcc of freedom dynamical 
systems. The continuous smooth evolution of frequency through the domain of initial conditions 
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does not exhibit the typical signature of a separatrix associated with an unstable periodic orbit. 
The range of the variation of the period of Atu is between 6,500 and 10,000 years: the shortest 
period occurs when the system is near Mode I and the longest period occurs near Mode II. The 
initial co-ordinate of the v Andromedae system (planets C and D) is indicated by the dashed 
vertical line; the oscillation period of Aw in this case is ^ 7, 100 years. 

Figures 6 and 7 are analogous to Figures 4 and 5, respectively, except they were obtained for 
the lower value of the normalized secular energy, iJgec = —0.4295. For this lower energy, we find 
that a qualitatively new feature appears in the phase space. In Fig. 6, we see the formation of 
a new libration center and one saddle-point in the domain of high eccentricities, where the effect 
of nonlinearity of the secular perturbations is strong. This feature can be associated with the 
advent of a nonlinear secular resonance zone, which is bounded by a zero-frequency separatrix 
and which occupies the domain about the Mode I fixed point at large eccentricities. This feature 
can be identified in the spectral map shown in Fig. 7, where the Aoj-frequency tends to at 
ei cos Aw = 0.36. 



4. Dependence on the Planetary Masses and Semimajor Axes Ratios 

In the previous sections, we have analyzed the secular behavior for one set of values the 
planetary mass and semimajor axis ratios: ai/a2 = 0.332 and mi/m2 = 0.461, values which are 
pertinent to the current v Andromedae planetary system. In this section, we will vary these two 
parameters and see the effect on the secular phase space structure. To see this dependence, we will 
plot the locations of periodic orbits and the domains of circulation or oscillation of Aw in the (61,62) 
representative plane of initial conditions. We will first consider the case when the outer planet is 
equal or more massive than inner one, that is m\/m2 < 1. Fixing the inner mass at mi = 1.83Mjup, 
we take three values of the outer mass: mi/m2 = 1/1, mi/m2 = 1/2 and mi/m2 = 1/4. Also, 
fixing the inner semimajor axis at oi = 0.83 AU, we take three values of the outer semi-major axis: 
02 = 1.5 AU, 02 = 2.5 AU and 02 = 5AU. 

The results are shown in Fig. 8, where the rows present the variation with mass parameter, 
for a fixed value of 01/02, while the columns present the variation with axis parameter, for a fixed 
value of mi/m2- Since the most significant changes occur at Aw = 0, only positive half-planes are 
presented in Fig. 8. In each plane, the boundaries of the domain of non-crossing planetary orbits are 
shown by thick inclined lines. The families of stable periodic orbits are shown by continuous black 
curves, whereas the boundaries of domains of Atu-oscillation are shown by dashed curves. Finally, 
the domains of the secular resonance are shown as shaded regions on the (ei,62)-planes. Note that 
the parameters (planet mass and semimajor axis ratios) used in construction of the mid-row and 
mid-column are very close to those of the current v Andromedae system shown in Fig. 3. 

Figure 8 gives a global view of all possible regimes of motion with respect to the secular 
angle Aw. When the planetary masses are equal and the mutual distances are small, almost the 
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Fig. 6. — Same as in Figure 4, but for iJgec = —0.4295. The erratic scatter of the points at 
eccentricities close to 1 is due to the loss of numerical accuracy. 
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Fig. 7. — Same as in Figure 5, but for -ffgec = —0.4295. 



whole phase space of the system is oscillations of the angle Aro around 0. In this case, as can 
be easily inferred from the geometrical pictures in Fig. 2, the oscillations of Aro around 180° will 
cover almost the whole negative half-plane of initial conditions (not shown here). The circulation 
occurs only in the domains of small eccentricities. With increasing mutual distance of the planets, 
the regions of circulation are also increasing. The same effect is observed when the mass of the 
outer planet increases. The secular resonance seems to be a common feature of the planetary 
system, particularly, when the masses and axes ratios are not close to 1. Its location on the (61,62) 
representative plane of initial conditions is always in the regions of the very-high eccentricities of 
the more massive outer planet. The orbit of the inner planet in this case can be nearly circular. 

Now we will consider the case when mass of the inner planet is larger that the mass of the 
outer planet, that is m\/m2 > 1. The geometrical picture of the (61,62) representative plane of 
initial conditions corresponding to this case is shown in Fig. 9. This plane was constructed with 
the same axes ratio used in Fig. 3, but the planetary masses were changed to nii = 3.97Mjup and 
777,2 = l-83Mjup. Since the axes ratio has not been changed, the energy level shown by the thick 
curve is the same level, Hgec = —0.41391, shown in Fig. 2. However, the structure of the ii'2-levels 
is different when compared to those in Fig. 3 because of the mass dependence of the variable K2 . 
The main geometrical difference is the location of the points B and C on the (ei,e2)-plane: now 
the point B is located in the negative half-plane and the point C in the positive half-plane. This 
means that the oscillation of Atu around occurs along the energy segment CD+ and around 180° 
along the energy segment AB. Finally, the initial conditions chosen in the intervals CA and BD~ 
correspond to circulating orbits with respect to Aw. 

The dependence of the secular behavior on the semimajor axes ratio in the case of mi/ 1712 = 2 
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Fig. 8. — Periodic orbits and domains of oscillation and circulation of Aro in the representative 
(ei,e2)-planes, for a range of planetary masses and semimajor axes. Rows: a fixed value of the inner 
semi-major axis oi = 0.83 AU and three values of the outer semi-major axis: 02 = 1-5 AU (top), 
a2 = 2.5 AU (middle) and 02 = 5 AU (bottom); Columns: a fixed value of the inner mass mi = 
l.SSMjup and three values of the outer mass, such as: mi/m2 = l(left), m\/m2 = 1/2 (middle) 
and mi/m2 = 1/4 (right). The shaded regions are domains of the nonlinear secular resonance. 



is shown in Fig. 10. This figure is similar to Fig. 8: it was constructed using the same three values of 
the outer semi-major axis 02 = 1.5 AU, 02 = 2.5 AU and 02 = 5AU, whereas the inner semimajor 
axis was fixed at 04 = 0.83 AU. As in Figure 8, the boundaries of the domain of noncrossing 
planetary orbits are shown by the thick inclined lines on each panel. The families of stable periodic 
orbits are shown by the thin continuous curve, whereas the boundaries of domains of Ac:7-oscillation 
are shown by the dashed curves. The main difference in the secular motion is that in this case the 
periodic orbits have lower values of 62 compared to the case with mi/m2 = 1/2, and the nonlinear 
secular resonance zone does not occur for all considered values of the axes ratio. Interestingly, the 
location of Mode I in this case is remarkably different from that predicted by linear secular theory. 
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Fig. 9. — Same as in Figure 2, but calculated with mi = 3.97Mjup and m2 = 1.83Mjup. 

5. Summary and Discussion 

The nonlinear secular perturbation analysis of the planar three body (two planet) problem pre- 
sented here has been motivated by the recent discovery of extra-solar planetary systems which have 
large orbital eccentricities and, in some cases, exhibit secular oscillations of their apsidal lines. The 
classical Laplace-Lagrange linear secular perturbation theory is insufficient for the precise analysis 
of such systems. We developed a semi-analytical approach, consisting of a numerical averaging over 
the short-period perturbations in the mutual interaction of the two planets to obtain the secular 
Hamiltonian. The Hamilton- Jacobi approach wc use shows that the secular Hamiltonian is sep- 
arable. (It is reducible to a one-degree-of-freedom integrable dynamical system for the canonical 
coordinate Aw.) We obtain the periodic orbits of this system, and we present a geometrical picture 
of the secular phase space of the two-planet system in terms of physical variables, the eccentricities 
and periastron longitudes. This phase space structure depends upon two parameters: the ratios 
of the planets' masses and their semimajor axes. The description of the secular variations of the 
system is completed by the initial conditions: the eccentricities and the relative longitude of the 
periastrons. Owing to the symmetries in the secular Hamiltonian, the phase space structure can be 
visualized in a two dimensional representative plane of initial conditions, 62 and ei cos Aw, where 
initial Alo is either or 180°. 

The locations of periodic orbits and the domains of different quasiperiodic motions (correspond- 
ing to the stationary solutions and oscillations or circulations of Aw) are then easily obtained in 
this plane of initial conditions, for given values of the parameters (mi /m2, 01/02). As expected from 
the linear secular perturbation theory for two planets, there exist two families of stable periodic 
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Fig. 10. — Same as in Figure 8, but for mi = 3.97Mjup and m2 = 1.83Mjup. The shaded regions 
are domains of oscillations of Aw around or 180°. 

orbits which describe the aligned and anti-aligncd orientations of the two planets' apsidal lines. 
Departures from the linear secular theory are readily visualized in this plane of initial conditions. 
In the domain of large eccentricities, we find a new feature which consists of a family of unstable 
periodic orbits. In this regime, the secular phase space contains a nonlinear resonance zone bounded 
by a zero-frequency separatrix. We emphasize that all of these results are obtained by means of a 
geometrical analysis based on the global constants of motion, and without time-expensive numerical 
integrations of the equations of motion. 

We applied this analysis to the specific case of the system of two outer planets, C and D, of v 
Andromedae, assuming edge-on co-planar orbits. We find that the nonlinear effects are significant 
in this case for eccentricities exceeding ~ 0.3 in the neighborhood of the aligned mode, and ~ 0.1 in 
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the neighborhood of the anti-ahgned mode. As noted aheady in previous studies of the dynamics 
of this system, the present best-fit initial conditions place this system in the domain of oscillations 
about the periodic orbit describing the aligned mode (Mode I, Atj7 = 0). Our analysis shows 
that the secular behavior of the v Andromedae system is stable over a large domain of initial 
eccentricities. On the representative plane of initial conditions, the domain of stable motion has 
an upper bound at the collision curve defined by close approaches between two planets; below the 
collision curve, a secular instability occurs in the vicinity of the nonlinear secular resonance located 
near ei = 0.6 and Aw = (Fig. 3). 

We have presented an analysis of the secular dynamics of the three-body planetary system 
over a wide range of the planetary mass and semimajor axis ratios. The systems always exhibit two 
main regimes of motion, characterized by circulation of Aw or its oscillation around or 180°. The 
existence of a third regime, the nonlinear secular resonance, depends on the values of the physical 
parameters. For example, when the planetary masses are equal and their mutual distances are small, 
the nonlinear secular resonance does not occur; however, under these conditions, the probability of 
finding the system in the oscillation (around or 180°) regime of motion is very close to 1. When 
both mass and semimajor axis ratios are far from 1, the domains of oscillation of Aw decrease, 
and the nonlinear secular resonance zone appears. The secular resonance feature discovered in our 
nonlinear analysis is presently only a theoretical curiosity, as no real system is known which can 
be associated with this dynamical state. However, the regime of planetary parameters where this 
feature is significant is not especially extreme, and it would not be surprising if such systems are 
identified in the future. 

The origin of the high planetary eccentricities of the v Andromedae planets C and D, as 
well as of many other exo-planetary systems, are an outstanding problem of high current interest. 
Two basic physical mechanisms that have been proposed to account for the high eccentricities 
invoke dynamical planct-planct interactions (Ford et al. 2001, Marzari & Weidenschilling 2002) 
or planet-disk interactions (Goldreich & Sari 2003); more complex scenarios involving planet-disk 
interactions moderated by resonant interactions between two planets have also been suggested 
(Chiang et al. 2002, Lee & Peale 2002, Malhotra 2002, Chiang & Murray 2002). It has been pointed 
out that in some multiple-planet systems the occurrence of apsidal oscillation, and especially the 
amplitude of that oscillation, may serve to constrain the physical mechanisms (Malhotra 2002). 
Unfortunately, the observational uncertainties of the planetary parameters of even the best-known 
exo-planetary system, v Andromedae, are still too large to define sufficiently precisely the amplitude 
of oscillation of Aw, and thus to favor one scenario over another. We await future improvements 
in the accuracy of the parameters of this system and others to more precisely determine their 
dynamics and thereby constrain their dynamical history. 

The method developed here is a non-linear generalization of the classical secular perturbation 
theory for the planar three body problem. In future work we hope to extend this to non-coplanar 

systems, and to apply the analysis to the secular dynamics of other planetary systems. We note 
that the Hamilton-Jacobi approach used here also provides a powerful means to explore the effects 
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of adiabatic evolution of the energy and angular momentum due to dissipative interactions with 
the nascent protoplanctary disk, identify adiabatic invariants, and thus possibly obtain additional 
insights into the dynamical history of some exo-planetary systems. 
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